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New methodology for verification of computational methods using unstructured grids is presented. The 
discretization order properties are studied in computational windows, easily constructed within a collection of 
grids or a single grid. The windows can be adjusted to isolate the interior discretization, the boundary dis- 
cretization, or singularities. A major component of the methodology is the downscaling test, introduced previ- 
ously for studying the convergence rates of truncation and discretization errors of finite-volume discretization 
schemes on general unstructured grids. Demonstrations of the method are shown, including a comparative ac- 
curacy assessment of commonly-used schemes on general mixed grids and the identification of local accuracy 
deterioration at intersections of tangency and inflow/outflow boundaries. Recommendations for the use of the 
methodology in large-scale computational simulations are given. 


I. Introduction 

There is an increasing reliance on computational simulations in aircraft design practices, supplementing traditional 
analytic and experimental approaches. Verification and validation methodologies 1 are being developed to ensure the 
correct application of these simulations. Verification methodologies for structured grids 2 are relatively well-developed 
in comparison to unstructured grids, especially grids containing mixed elements or grids derived through agglomera- 
tion techniques. The summary of the latest of three Drag Prediction Workshops 3 illustrates the problems associated 
with assessing errors in practical complex-geometry/complex-physics applications. Current practices tend to compare 
relative errors between methods and experimental results rather than absolute errors. The motivation for this paper 
was to propose verification methodologies to understand precisely the expected code performance in such large-scale 
computational endeavors. 

In an earlier paper, 4 analytical and computational methods for evaluating the accuracy of finite-volume discretiza- 
tion (FVD) schemes defined on general unstructured meshes were presented. The study corrected a misconception 
that the discrete solution (discretization) accuracy of FVD schemes on irregular grids is directly linked to convergence 
of residuals evaluated with the exact solution (truncation errors). The definitions and properties of discretization and 
truncation errors are discussed subsequently in this paper and also elsewhere. 4 

The main computational tool introduced is a downscaling (DS) test. Performed for a known exact or manufactured 
solution, the test consists of a series of inexpensive computational experiments that account for local properties of 
the studied scheme; it provides estimates for the convergence orders of the discretization and truncation errors by 
comparing errors obtained on different scales without any restrictions on the grid structure. It can be tailored to 
study the discretization accuracy in the interior, at the boundary, and/or in vicinity of singularities. Analysis methods 
predicting the performance of DS tests have also been developed. 

The overall testing process is shown in Figure 1, summarizing relations between convergence orders of discretiza- 
tion and truncation errors predicted by these analytical estimates, observed in DS tests, and observed in global grid- 
refinement computations. The analytical estimates predict error convergence observed in DS tests; the estimates are 
conservative because they do not account for possible error cancellation occurring on regular (i.e., mapped) grids. The 
DS tests provide exact truncation error convergence, but can be optimisitic predictors of discretization error because 
they do not account for possible global discretization-error accumulation. 
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The entries in Figure 1 are arranged from lowest to highest compu- 
tational cost. Unfortunately, the less expensive estimates are are more 
difficult to interpret correctly. For example, grid refinement in large- 
scale 3D computations is quite expensive. But it is quite simple to as- 
certain attainment of design order in grid refinement if an exact solution 
is available. In monitoring truncation errors, the solutions need not be 
determined, only residuals need to be evaluated with the manufactured 
solution. Truncation error assessment is inexpensive but is more diffi- 
cult to interpret; when interpreted correctly, it is a powerful verification 
tool. Because DS-test estimates of discretization errors neglect certain 
error accumulation mechanisms, the extent to which they can be trusted 
as predictors of global accuracy is an area of current research. In any 
case, the DS-tests are always optimistic predictors of discretization er- 
ror, so if they fail, there is certainly an error in the formulation. Thus, 
the hierarchy of analytical and computational estimates can be used to 
complement current verification practices in large-scale simulations. 

In large-scale grid-refinement studies, extensive amounts of data are 
involved and integral norms often do not provide sufficient information 
to isolate the source of errors. An attractive feature of the approach advo- 
cated here is that tests can be tailored to pinpoint regions/solutions/grids 
of interest. Also, in DS-testing, very small meshes sizes can be used to ensure that testing is within the asymptotic 
convergence range (where the leading-order terms dominate). 

The purpose of this paper is to extend the verification methodology towards the practical verification of large-scale 
general unstructured grid formulations. The overall approach is applied through windows defined on a computational 
domain; grid-refining and/or DS tests are performed in each window. Demonstrations are shown, including a compar- 
ative accuracy assessment of commonly-used FVD schemes on general unstructured grids of mixed type. The local 
accuracy deterioration at intersections of tangency and inflow/outflow boundaries is identified using tailored DS tests. 

II. Elements of Analysis 

The FVD schemes are derived from the integral form of a conservation law 

f (F • n) dT = f j (f - S) (El, (1) 

r n 

where f is a forcing function independent of the solution. S' is a solution-dependent source function, O is a control 
volume with boundary F, n is the outward unit normal vector, and F is the flux vector. The main accuracy measure of 
any FVD scheme is the discretization error , E d , defined as the difference between the exact continuous solution, Q, 
to the differential conservation law 

VF = / - S (2) 

and the exact discrete solution, Q h , of the discretized equations (1) 

E d = Q~ Q h . (3) 

A common approach to evaluate the accuracy of discrete schemes is to monitor the convergence of truncation 
errors. Truncation error, E t , measures the accuracy of the discrete approximation to the differential equations (2). 
It is defined as the difference between the differential equations and their discrete representation 5 6 and is found by 
computing the normalized discrete residuals after substituting the exact solution for the discrete solution. For schemes 
approximating integral conservation laws, the truncation error is usually defined from a time-dependent standpoint. 7 - 8 
In the steady-state limit, after substituting the exact solution Q into the normalized discrete equations (1), the truncation 
error is defined as 


Analysis Estimates 

• conservative for truncation errors 

• conservative for discretization errors 

DS-Test Estimates 

• exact for truncation errors 

• optimistic for discretization errors 

Grid-Refinement Computations 


Figure 1. Overall testing process. 
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where F h is a reconstruction of the flux F at the boundary F, Cl is the measure of the control volume, 

\n\ = jj dto, (5) 

n 

f h and S h are, respectively, approximations of the forcing function f and the source function S on Cl, and the integrals 
are computed according to some quadrature formulas. 

Convergence of truncation errors as a FVD accuracy measure on structured (regular) grids is well justified; the 
truncation errors converge as 0(h p ) on sequences of refined meshes, where h is a characteristic mesh size and p is 
the design discretization-accuracy order of the method. However, the truncation-error convergence is often misleading 
for unstructured grids. Several studies 9,10 noted that 2 nfi - order convergence of truncation errors for some commonly 
used FVD schemes can be achieved only on grids with a certain degree of geometric regularity. Examples published 
elsewhere 3 4 and in this paper show that the design-order discretization-error convergence can be achieved even when 
truncation errors exhibit a lower order of convergence or, in some cases, do not converge at all. For some formally 
inconsistent FVD schemes (truncation errors do not converge), it can be rigorously proved that the discretization errors, 
in fact, converge. 1 1 One interpretation of these results is that the FVD schemes on truly unstructured grids approximate 
a perturbed Equation 2 with an additional oscillatory forcing function. The mean value of the perturbation computed 
over an arbitrary portion of the domain still approaches zero with the design order. 

The manufactured solutions used herein are of two types, either simple analytic functions (collections of poly- 
nomials or sines) or exact solutions. The corresponding forcing functions are found by substituting these solutions 
into the continuous governing equations and boundary conditions. The intent of the approach is to facilitate testing of 
discretizations and boundary conditions in situ for large-scale computations; this is possible with slight modifications 
of most boundary conditions, e.g., evaluating no-slip conditions with a specified wall velocity instead of the usual zero 
velocity condition. Likewise, in the farfield, the exterior conditions are taken from the exact solution rather than the 
usual assumption of constant freestream exterior conditions. Not all boundary conditions are amenable to such a mod- 
ification (e.g., inviscid tangency) and for these we use exact (or manufactured) solutions associated with a particular 
geometry. An alternative is the mapping construction used by Bond et al. 12 


III. Consistent Refinement 




The general FVD approach requires par- 
titioning the domain into a set of non- 
overlapping control volumes and numeri- 
cally implementing equation (1) over each 
control volume. Two types of FVD schemes 
are considered: node-centered schemes, in 
which solution values are defined at the 
mesh nodes, and cell-centered schemes, in 
which solutions are defined at the centroids 
of the control volumes. In the 2D examples 
considered here, the primal meshes are com- 
posed of triangular and quadrilateral cells; 
in 3D computations the cells are tetrahedral, 
prismatic, or hexahedral. The median-dual F| S ure 2 - 
partition 13, 14 used to generate control vol- 
umes for the node-centered discretization is 

illustrated in Figure 2(a) for two dimensions. These non-overlapping control volumes cover the entire computational 
domain and compose a mesh that is dual to the primal mesh. For cell-centered FVD schemes, the primal cells serve as 
control volumes (Figure 2(b)). 

The discrete solution is represented as a piecewise linear function defined either within primal or dual cells. The 
discretizations are applied at a sequence of refined grids satisfying the consistent refinement property. This property 
requires the characteristic distance across primal and dual cells to decrease consistently with increase of the total 
number of degrees of freedom, N. The characteristic distance should tend to zero as N~ 1 where d is the number 
of spatial dimensions. The property enables a meaningful assessment of the asymptotic order of error convergence. 
In particular, on 3D unstructured meshes satisfying the consistent refinement property, the discretization errors of 


(a) Node-centered. 


Control volumes (shaded area) for FVD schemes; grid points are solid 
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second-order FVD schemes are expected to be proportional to N ~ 2 / 3 . 

An equivalent mesh size based on the degrees of freedom is defined as Tin = N~ 1 ^ d . An equivalent mesh size 
based on a characteristic distance is defined in terms of norms of the local control-volume function, i.e., hy = \ V 1 / ' , | 
where 1 1 • 1 1 is a norm of choice. For a consistently refined mesh, hy is a linear function of h ..y for any physical regions. 
The assessment of consistent refinement is purely geometric and could be done automatically by inspecting the mesh 
over local subsets of the domain. The procedure is envisioned to be most useful during the grid generation phase to 
identify and and repair regions where the grids are not consistently refined. 



(c) Coarse grid (nearfield). (d) Fine grid (nearfield). 

Figure 3. Partial view of surface grids on symmetry plane and sphere. 

To illustrate the concept, we analyze three unstructured tetrahedral grids generated around a sphere; the grids are 
composed of 25,473 nodes, 82,290 nodes, and 328,463 nodes. In Figure 3, farfield and nearfield views of the coarsest 
and finest surface grids are shown. In Figure 4, variations of hy based on the L \ - and L <*, -norms of V 1 / 3 are shown 
versus the equivalent mesh size Tin = 1//3 , each normalized by the value on the coarsest grid. A consistently refined 

mesh variation is denoted by a dashed line in the figure. Based on the Ly norm, hy is linear but the hy computed 
with the Loo-norm shows that the mesh is not consistently refined. Examination of the grids in Figure 3 confirms 
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Equivalent Mesh Size (Degrees of Freedom) 


Figure 4. Consistent refinement check using normalized equivalent mesh sizes. 



W . I W . I 

Equivalent Mesh Size (Control Volumes) 


Figure 5. Variation of Li-norm of error in pressure with grid refinement. 


that the mesh near the farfield boundary is not consistently refined. Inviscid incompressible flow equations discretized 
with a second-order node-centered FVD scheme have been solved on these grids. The Li-norms of the errors in 
pressure, shown versus hy in Figure 5, converge with second order, probably because the solution variations are much 
larger near the surface than near the farfield boundary. With continued inconsistent refinement, the discretization error 
convergence would eventually degrade. 

The discretization error on the coarsest grid using a cell-centered second-order FVD scheme is also shown in 
the figure. The cell-centered scheme has approximately six times as many degrees of freedom as the node-centered 
method on the same grid. The discretization error is comparable for the two schemes when the number of degrees 
of freedom are matched. Other results (not shown) generated with consistently-refined meshes over segments of the 
sphere confirm that the discretization accuracy of cell-centered and node-centered grids is comparable for comparable 
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hy. Note this result differs from the transonic results reported by Levy and Thacker 15 and Mavriplis 16 in which the 
errors between node-centered and cell-centered results are correlated by matching the number of degrees of freedom 
on the surface. 

IV. Downscaling (DS) test 

The DS test is designed to predict the local convergence orders of discretization and truncation errors observed 
in large-scale computations. The test provides exact estimates for truncation error convergence, but does not account 
for possible discretization-error accumulation. On truly unstructured multidimensional grids, the local discretization 
errors are dominant; therefore, the DS-test estimates of the discretization-error convergence orders are expected to be 
accurate as well. 

The DS test employs numerical computations on a sequence of contracted domains zooming toward a focal point 
within the original computational domain (Figure 6). The choice of the manufactured solution and the focal-point 
location can be tailored to study accuracy of a typical interior discretization, a boundary discretization, or a specific 
solution singularity. Because the number of degrees of freedom in the DS-test domains is held (approximately) fixed, 
one can study convergence properties on grids with characteristic mesh sizes much smaller than those affordable for 
practical computations with a globally-refined grid sequence. 

There are at least two pos- 
sible strategies for grid gen- 
eration on these contracted 
domains. The first strategy 
is termed “scaled grid” (Fig- 
ure 6(a)). With this strategy, 
the first (coarsest) computa- 
tional domain is debited as 
a subdomain of the investi- 
gated global mesh containing 
the focal point; other (bner) 
domains and their mesh pat- 
terns are derived by scal- 
ing down this brst domain 
(e.g., repeatedly multiplying 
all the distances from the fo- 
cal point by a given factor, 
say, 1/2 or 1/4). An in- 
dependent grid (Figure 6(b)) 
can be generated on each 
domain, assuming that the 
consistent rehnement prop- 
erty is satisbed, i.e., the char- 
acteristic distance across a f a ) Scaled § rid - ( b ) Independent grid generation (ac- 

grid cell is scaled down with ^ v : 

the same rate as the diam- 
eter of the contracted do- 
mains The scaled-orid ap- Figure 6. Illustration of DS computational domains. Black bullets mark the focal points; white 
. . ... squares mark the interface between the interior and the DS-test domain, 

proach is especially useful tor 

studying interior discretiza- 
tions and straight boundaries. It is impractical for studies near a general (discretely debited) curvilinear boundary 
because the physical boundary shape should be preserved on each grid in the DS sequence. 

In the DS test, the studied FVD scheme is supplemented with a set of boundary conditions at the interface between 
the interior and the DS domain (marked with white squares in Figure 6); overspecihcation from the known manufac- 
tured solution is a typical choice. In studying accuracy of discretized boundary conditions, i.e., when the focal point is 
at the physical boundary, the physical conditions are implemented at the boundary surface; overspecihcation can still 
be applied at the interface between the interior and the DS domain (see sketches in Figure 6). The freedom to choose 
the manufactured solution, the focal point, the shape of domains, and the type of interface boundary conditions greatly 
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simplifies DS testing. 

The DS test evaluates local truncation and discretization-error convergence orders by comparing errors obtained 
in computations on different scales. The convergence of truncation errors (in the Loo-norm) observed in global grid- 
refinement computations will be bounded by the worst DS-test estimate, provided DS-testing samples all representative 
regions. Global convergence in integral norms, e.g., L \ -norm, may be higher order because these norms are less 
sensitive to fluctuations occurring locally. 

On grids with an increased regularity, the DS test may overestimate convergence of the discretization errors. Some 
global error accumulation can occur on such grids; because of its local nature and overspecified interface boundary 
conditions, the DS test is incapable of accounting for the global error accumulation. Our experience shows that 
on truly unstructured multidimensional grids, errors produced locally dominate globally accumulated errors, and the 
discretization-error estimates obtained in a DS test become accurate. 


V. Windowing 

In order to provide a framework for assessing performance of practical large-scale codes for specific computations, 
we introduce the concept of windowing. A window is an arbitrarily-shaped subdomain within the computational 
domain (see sketch in Figure 7) serving as a reference frame for testing. Each test captures an entry from the three 
features affecting convergence of local errors: (1) problem-related features; (2) solution-related features; and (3) 
discretization/grid-related features. 

The problem-related fea- 
tures are determined by the 
scope of required compu- 
tations. Specifically, the 
features include the interior 
governing equations, vari- 
ous types of boundary con- 
ditions (e.g., inflow, outflow, 
tangency, no-slip, symme- 
try), and the geometrical fea- 
tures characterizing bound- 
aries (e.g., smooth flat bound- 
ary, curved boundary, sharp 
corners). To address the 
problem-related features, the 
windows should be placed in 
representative locations (in- 
terior, boundaries, corners, 
etc.) 

The solution-related fea- 
tures account for variations 
in the solutions typically en- 
countered, including smooth 
flows, shocks, stagnation re- 
gions, vortices, boundary lay- 
ers, recirculating flows, etc. Each feature should be represented by a specific choice of the manufactured solution. 

The discretization/grid-related features concern variations in meshes and discretization schemes. The features in- 
clude the interior discretization scheme, discretization of boundary conditions, grid composition (e.g., combinations of 
advanced-layer (prismatic) regions with interior tetrahedral regions), approximation of geometry (flat panel or higher- 
order approximation), etc. Interfaces between regions with different meshes represent separate grid-related features. 
Allowed grid singularities, such as hanging nodes, degenerated cells, and similar ones should also be considered as 
separate grid-related features. 

To verify a large-scale code for certain computations, each representative triplet of features requires a designated 
test; convergence of truncation and discretization errors observed in all representative tests has to be understood and 
accepted as satisfactory. For example, in the interior of structured grids, the truncation errors of FVD schemes for 
smooth manufactured solutions are expected to satisfy the design order. On unstructured grids, the norms of trun- 



Figure 7. Sketch of possible windows superimposed on an unstructured grid. Regions denoted by 
solid line are grid-scaled windows (with black focal points); regions denoted by dashed line are win- 
dows preserving body geometry (with gray focal points). 
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cation errors often converge slower than corresponding norms of discretization errors; the truncation error conver- 
gence is allowed to be 1 order lower than the design order for inviscid-dominated equations and 2 orders lower for 
viscous-dominated terms. 4 The average truncation errors observed in grid-refinement computations performed within 
a window are expected to converge with the design order. 

Establishing the discretization-error convergence order in grid-refinement computations is not practical because 
discrete solutions must be computed on grids with prohibitively many degrees of freedom. The DS tests enable such 
calculations because the number of points on each grid is kept (approximately) constant. One should interpret the DS 
test results carefully; in particular, on structured grids, convergence of discretization errors observed in DS tests is 
expected to be higher-order than that observed in grid-refinement computations. On unstructured grids, where locally- 
generated errors are dominant, the discretization-error convergence observed in DS tests is of the same order as the 
convergence observed in global grid-refinement computations. In any case, if the observed convergence of either 
discretization or truncation errors is slower than expected, this is an unambiguous indication of deficiencies in either 
formulation or implementation. 

In some cases, these deficiencies may be found acceptable, for example, when large discretization errors are 
generated and remain local, without affecting integral norms of the errors computed over the entire domain. 

VI. Example 1: Two-dimensional Laplace equation 

To illustrate applications of DS tests, we, first, consider the two-dimensional Laplace equation, as a model of the 
diffusion terms in the Navier-Stokes equations, 

A U = /, (6) 

subject to Dirichlet boundary conditions. The equations are discretized with a second-order edge-based node-centered 
FVD scheme defined on a series of random mixed-element grids composed of triangles and quadrilaterals. Details of 
the discretization can be found in 417 The manufactured solution and forcing term are taken as U = [(sin(7ra:)) 2 + 
(sin(7rj/)) 2 ]/2, / = — 27r 2 [l — (cos(7ra;)) 2 — (cos(7n/)) 2 ]. 

For illustration purposes, 
the computations performed 
in DS test are compared with 

global grid-refinement com- 2 

putations. For global grid re- 
finement, each grid is formed 

from an underlying struc- 15 

tured quadrilateral grid (Fig- 
ure 8). In terms of a po- i 

lar, (r, 0), coordinate system, 
the grid extent is defined as 
9 G [7t/6,7t/ 3] in the circum- Figure 8. A typical mixed-element unstructured grid generated with random splitting and random 
ferential direction and r G perturbation of the underlying quadrilateral grid. 

[1, 2.2] in the radial direction. 

The decision to split (or not to split) each structured quadrangle into triangles is determined randomly; approximately 
half of the quadrilaterals are split. In addition, the interior grid points are perturbed from their original position by ran- 
dom shifts in the range [— %/ 2/6, x/2/6] of the local mesh size in the radial direction. The sequences of globally refined 
grids are generated with 2 n+3 + 1 points in both the radial and circumferential directions, where n = 0, 1, 2, 3, 4. 
The sequences of DS grids are generated from a grid with 17 points in both the nominal radial and circumferential 
directions and downscaled about the center of the domain by a factor 2 -s , where s = 0, 2, 4, 6, 8. The grid topology 
remains unchanged. 

The Li norms of truncation and discretization errors are shown in Figure 9 versus an equivalent mesh size pa- 
rameter, hy. Although not shown, error convergence rates in the norm are the same as the Li-norm rates. In 
grid-refinement computations, the truncation errors remain 0(1) and the discretization errors converge with 2" d -order, 
precisely predicted by the DS test. The reason for the 0(1) convergence of truncation errors is grid irregularity stem- 
ming from the usage of unstructured grids. The literature frequently associates 0(1) convergence of truncation errors 
on irregular grids as an indication of a inconsistent scheme that never converges to the exact result; 9,18 this example 
clearly shows that this is not necessarily the case. The convergence orders of both truncation and discretization errors 
are within the ranges allowed for second-order FVD schemes. 
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Equivalent Mesh Size 



Equivalent Mesh Size 


(a) DS test (b) Grid refinement computations 


Figure 9. Convergence of the discretization and truncation errors for the Laplace equation solved on irregular mixed-element unstructured 
grids. 


VII. Example 2: Two-dimensional incompressible Euler equations 

In this section we consider incompressible inviscid equations in the interior and next to the curved tangency 
boundary . Inviscid fluxes for conservation of mass and momentum are defined as 

/ 3u f3v 

F = fi + gj= u 2 +p i+ uv j, (7) 

uv v 2 + p 

where the vector of unknowns, Q = [it, v. p\, includes the Cartesian velocities and the pressure; f3 is an artificial 
compressibility parameter, 17 taken as /3 = 1 here. 

Two common FVD schemes with design second-order accuracy are investigated: a edge-reconstruction median- 
dual node-centered scheme and a cell-centered scheme. The node-centered FVD scheme uses the least square method 
for gradient reconstruction and integration over the control-volume boundaries employing fluxes evaluated at the edge 
medians; details of the discretization can be found elsewhere. 4, 17 The cell-centered FVD scheme also employs the least 
square method for gradient reconstruction. 13 Numerical tests presented are performed for 2D inviscid incompressible 
non-lifting flow around a cylinder of unit radius centered at the origin. The analytical solution for this problem is well 
known. 4 

The first set of tests is performed to study accuracy of the interior discretization. The computational domain is 
shifted away from the surface of the cylinder: 1.5 < r < 4, 27r/3 < 0 < Itt/ 3. The two FVD schemes are studied 
on random triangular and random mixed-element grids. Examples of unstructured grids derived from an underlying 
structured grid are shown in Figure 12. Grid randomization is introduced through random splitting (or not splitting) of 
structured quadrilateral cells. Each cell has equal probabilities to introduce either of the two diagonal choices or, for 
mixed-element grids, no diagonals. 

For each formulation, grid-refinement and DS tests are performed. In global grid refinement computations, the 
underlying structured grid is refined by doubling the number of intervals in the radial and angular directions. Ran- 
domization is introduced independently on each scale. The inflow boundary conditions are enforced at the boundary 
corresponding to the external radius; outflow conditions are enforced at all other boundaries. In the DS test, the coars- 
est 9 x 9 grid is scaled down around the point r = 2.75, 6 = tt by multiplying all angular and radial differences from 
this point by a factor of 0.5. Table 1 summarizes the convergence of discretization and truncation errors observed in 
these tests. The convergence orders are the same in all norms and for all variables and equations. The results con- 
firm capabilities of the DS test to provide accurate estimates for error convergence in grid-refinement computations. 
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The observed convergence rates indicate that the edge-based node-centered FVD scheme is second-order accurate on 
triangular grids, but only first-order accurate on mixed-element grids; the cell-centered formulation is second-order 
accurate on all studied grids. 


Formulation 

DS test 

Grid-refinement computations 

Discr. Error 

Trunc. Error 

Discr. Error 

Trunc. Error 

Node-centered, 
random triangular grid 

0(h 2 ) 

O(h) 

0(h 2 ) 

Oih) 

Node-centered, 
mixed-element grid 

O(h) 

0(1) 

0{h) 

0(1) 

Cell-centered, 
random triangular grid 

o(h 2 ) 

O(h) 

o(h 2 ) 

0(h) 

Cell-centered, 
mixed-element grid 

o(h 2 ) 

O(h) 

o(h 2 ) 

0(h) 


Table 1. Convergence of discretization and truncation errors for various unstructured grid formulation of the 2D inviscid incompressible 
equations on an inflow/outflow computational domain. 


There are many ways to recover second-order accu- 
racy with the node-centered FVD scheme. For example, 
second- and third-order node-centered schemes have been 
demonstrated with face-reconstruction techniques for flux 
evaluation. 4 

Although not shown, we have implemented a central 
version of the edge-reconstruction node-centered scheme 
and tested it for various unstructured grids. First-order 
convergence of discretization errors has been observed on 
mixed-element and random quadrilateral grids. The results 
are consistent with some previously published results 19 
and contradict others. 9 In the latter reference, 0(1) con- 
vergence of discretization errors on randomly-perturbed 
quadrilateral grids with a central scheme was reported; 
in-depth investigation of the discrepancies has been re- 
ported. 4 

Another series of tests has been performed to study ac- 
curacy of the FVD schemes at the curved tangency bound- 
ary; both schemes use isotropic triangular grids approxi- 
mating the curved tangency boundary by straight segments 
linking grid nodes located at the physical boundary. The 
approximation is illustrated in Figure 13 (a). The discrete 
tangency is enforced weakly over the straight segments. 

A sequence of random triangular grids is generated at 
the top of the cylinder (1 < r < 2.2, 7t/3 < 9 < 27t/ 3); 
a grid example is shown in Figure 13 (b). Figure 14 il- 
lustrates convergence of the Li norm of truncation and 
discretization errors in DS tests performed with the node- 
centered edge-reconstruction FVD scheme. The left figure 
(a) exhibits convergence observed in the DS test with the 
focal point in the middle of the tangency boundary; the 
right figure (b) shows results for the DS test with the focal 
point next to the inflow/tangency corner. See sketches in 
Figure 14. 



Figure 10. Random triangular. 



Figure 11. Random mixed. 

Figure 12. Typical triangular and mixed-element unstructured 
grids defined on a computational domain shifted away from the sur- 
face of the cylinder. 


The 2 " ' , -ordcr convergence of discretization errors and the l st -order convergence of truncation errors demonstrated 
in the interior-tangency DS test comply with expectations for a second order accurate scheme. Convergence deterio- 
ration is clearly observed in the DS test performed with the inflow/tangency boundary conditions, indicating local loss 
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(a) Straight-segment approxima- (b) Random triangular grid around the top of the cylinder, 

tions to curved tangency bound- 
ary (dashed line). 


Figure 13. Boundary approximation and grids for DS test of local boundary conditions. 




(a) DS test: interior tangency boundary condition. (b) DS test: inflow/tangency boundary conditions. 


Figure 14. Convergence of the L\ norm of rr— momentum truncation and discretization errors observed in DS tests performed on random 
triangular grids surrounding the top tangency boundary of the unit cylinder; Dashed and dashed-dot lines denote 1 st and 2 nd order error 
variations; open squares in sketch denote boundaries with overspecification. 


of second-order accuracy. This local accuracy deterioration is explained and repaired elsewhere. 4 The convergence of 
the L i norms of the discretization errors in the corresponding grid-refinement test show the second-order convergence, 
while the norms of the errors converge with the first order. These tests can serve as examples that local accuracy 
deterioration can be acceptable, if the cause and effect on discretization errors are fully understood. 

Analogous DS tests performed for the 2D cell-centered FVD scheme showed second-order convergence of dis- 
cretization errors at the interior tangency and at the inflow/tangency corner. 

VIII. Recommendations on verification procedure 

In this section, we provide recommendations on choosing relevant tests to verify a large-scale code for a certain 
type of computation; the illustrative examples are motivated by the recent Drag Prediction Workshops. 3 

There are two preliminary tests concerned with truncation-error computations (no need to compute discrete solu- 
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tions), which are useful for confirming consistency of the investigated FVD scheme. The first test is performed for a 
smooth manufactured solution on structured (e.g., Cartesian) consistently-refined meshes; design-order convergence 
of truncation errors is expected. The second test is performed for a conservation law equation and a manufactured 
solution that produces linear fluxes, for example, mass conservation with constant density and linear velocities, or mo- 
mentum conservation with constant density and velocity and linear pressure. Second (or higher) order FVD schemes 
are expected to exhibit zero truncation errors for the equations associated with linear fluxes on any mesh. 

Assuming the FVD scheme passed the consistency tests, the first step toward forming a library of the DS tests 
for a large-scale code verification is to formulate a list (as complete as possible) of relevant problem-, solution- and 
discretization/grid-related features. The following list has been compiled for a mixed-element unstructured-grid solver 
considered for computations of a viscous flow around an airfoil. 

• Problem-related features include Navier-Stokes equations with a given set of parameters, such as Mach and 
Reynolds numbers; turbulence model; far-field, symmetry, and no-slip boundary conditions; straight or smoothly 
curved profiles for the far-field and symmetry boundaries, smooth and discontinuous boundary profiles for the 
airfoil surface. Each problem-related feature is addressed in DS testing by choosing an appropriate computa- 
tional window. 

• Solution-related features include smooth non-degenerated flow, stagnation flow, vortex, shock, boundary layer, 
and flow separation. Various solution features are allowed to interact. Each solution-related feature is addressed 
in DS testing by choosing an appropriate manufactured solution. 

• Discretization/grid-related features include the interior FVD scheme, boundary discretization scheme, advanced- 
layer prismatic meshes within the boundary layers, and general tetrahedral meshes in the exterior. Interfaces 
between the regions with different meshing and mesh singularities should be considered as separate grid-related 
features. Each discretization/grid-related feature is addressed in testing by constructing the grid (grid-refinement 
generally requires additional grid generation whereas DS-test may not) and by applying appropriate discrete 
equations. 

A designated DS test should be designed for each relevant triplet of features, one from each group. Not all triplets 
are relevant; for example, there is no need to test combination of far-field boundary and a boundary-layer solution. 
In many cases, formulation deficiencies can be found by observing slow convergence of truncation errors. For other 
cases, e.g., second-order unstructured FVD schemes for viscous equations, the truncation errors are allowed to 0(1) 
convergence; only discretization error convergence can then be decisive. 

As examples, let us consider the DS tests recommended for verifying the interior formulation for smooth non- 
degenerated solutions. A computational window is placed away from all physical boundaries and a smooth non- 
degenerated flow is chosen as the manufactured solution. For such tests, the requirements for convergence of truncation 
and discretization errors have already been formulated: truncation errors are expected to be bounded and discretization 
errors are expected to converge with second order. Note, that on inviscid scales, where convection dominates diffu- 
sion, the truncation errors are expected to converge with first order. At least four basic combinations of nonsingular 
meshes should be considered: (1) general prismatic meshes, (2) general tetrahedral meshes, (3) random mixed-element 
meshes, and (4) meshes with a smooth interface between the prismatic and tetrahedral regions. If certain mesh singu- 
larities (e.g., hanging nodes, zero-volume elements, other types of elements beside triangular prisms and tetrahedrons) 
are allowed, they should be tested in separate DS tests, usually in combination with the four basic nonsingular meshes. 

A similar set of DS tests is suggested to verify the interior formulation for a vortex solution. Note that, for inviscid 
equations at stagnation, the convergence of discretization errors of velocity components tends to degenerate by one 
order. 4 This degeneration may or may not be noticed depending on the flow Reynolds number. Even if observed, the 
increased discretization error may stay local and not affect convergence of the Li-norms of the discretization errors. 

For verifying the formulation for smooth solutions in the vicinity of a smooth surface, one has to place the focal 
point at the surface and run DS tests with general prismatic meshes and manufactured solutions representing boundary- 
layer flow, stagnation flow, and separated flow. For testing smooth solutions around sharply angled parts of the airfoil 
surface, the same manufactured solutions should be DS-tested on general mixed-element meshes. We have explored 
only a subset of the recommended practices to date. In particular, the expected asymptotic behavior for discontinuous 
solutions has yet to be addressed. 


12 of 13 


American Institute of Aeronautics and Astronautics 



IX. Concluding Remarks 


New concepts to verify large-scale computational codes using unstructured grids are presented. The method in- 
troduces computational windows, easily constructed within a collection of grids or a single grid, to examine order 
properties of the scheme over the entire domain or in regions isolating the interior discretization, the boundary dis- 
cretization, or singularities. A major component of the methodology is the downscaling (DS) test, introduced earlier for 
studying the convergence rates of truncation and discretization errors of finite-volume discretization (FVD) schemes 
on general unstructured grids. Demonstrations of the method are shown, including a comparative assessment of the 
accuracy from commonly-used FVD schemes, on general mixed grids. Local deterioration of accuracy at intersections 
of tangency and inflow/outflow boundaries are identified using tailored DS-tests. 
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Towards Verification of Unstructured-Grid Solvers 
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New methodology for verification of large-scale computational methods using unstructured grids is pre- 
sented. The discretization order properties are studied in computational windows, easily constructed within 
a collection of grids or a single grid. The windows can be adjusted to isolate the interior discretization, the 
boundary discretization, or singularities. A major component of the methodology is the downscaling (DS) test, 
introduced earlier for studying the convergence rates of truncation and discretization errors of finite-volume 
discretization (FVD) schemes on general unstructured grids. Demonstrations of the method are shown, in- 
cluding a comparative assessment of the accuracy of a common edge-reconstruction FVD scheme, known to be 
second-order accurate for inviscid equations on triangular grids, with alternative second and third order accu- 
rate face-reconstruction schemes on general mixed grids. Deterioration of the local accuracy at intersections 
of tangency and inflow/outflow boundaries are identified and repaired using DS tests. Additionally, analysis 
and computations are presented demonstrating that the usual first-order discretization of convection in the 
turbulence equations does not degrade the order of accuracy of the meanflow equations. 


I. Introduction 

There is an increasing reliance on computational simulations in aircraft design practices, supplementing traditional 
analytic and experimental approaches. Verification and validation methodologies 13 are being developed to ensure the 
correct applicability of these approaches in practical applications. Verification methodologies for structured grids 12 
are relatively well-developed in comparison to unstructured grids, especially grids containing mixed elements or grids 
derived through agglomeration techniques. The summary of the latest of three Drag Prediction Workshops 10 illustrates 
the problems associated with assessing errors in practical complex-geometry/complex-physics applications. The mo- 
tivation for this paper was to propose verification methodologies more applicable to such large-scale computational 
endeavors. 

In an earlier paper, 1 analytical and computational methods for evaluating the accuracy of finite-volume discretiza- 
tion (FVD) schemes defined on general unstructured meshes were presented. The study corrected a misconception that 
the discretization accuracy of FVD schemes on irregular grids is directly linked to convergence of truncation errors. 
In some cases, design-order discretization accuracy on general irregular unstructured grids was demonstrated even 
with no convergence of the truncation error. In fact, on truly unstructured grids, the convergence order of truncation 
errors is typically less than the discretization-error convergence order by 1 for inviscid equations and by 2 for viscous 
equations. Methods were introduced to monitor the truncation error convergence and assess if the design order of the 
discretization accuracy can be achieved. 

The main computational tool introduced for evaluating the accuracy of discretization schemes is a downscaling 
(DS) test. Performed for a known exact or manufactured solution, the test consists of a series of inexpensive compu- 
tational experiments that account for local properties of the studied scheme; it provides estimates for the convergence 
orders of the discretization and truncation errors by comparing errors obtained on different scales without any restric- 
tions on the grid structure. It can be tailored to study the discretization accuracy in the interior, at the boundary, and/or 
in vicinity of singularities. Analysis methods predicting the performance of DS tests have also been developed. 

The purpose of this paper is to apply the verification methodology towards the practical verification of large-scale 
general unstructured grid formulations. The overall approach is applied through windows defined on a computational 
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t National Institute of Aerospace (NIA), 100 Exploration Way, bdi skin @ nianet.org . This research was supported by the National Institute of 
Aerospace under NASA LaRC Cooperative Agreement NCC 1-02043 
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domain; grid-refining and/or DS tests are performed in each window. Demonstrations are shown, including a com- 
parative assessment of the accuracy of a common edge-reconstruction FVD scheme with alternative second and third 
order accurate face-reconstruction schemes on general unstructured grids of mixed type. Improvements are made to 
the local accuracy at intersections of tangency and inflow/outflow boundaries using insights from the DS tests. Addi- 
tionally, analysis and computations are presented demonstrating that the usual first-order discretization of convection 
in the turbulence equations does not degrade the second-order accuracy of the meanflow equations. 

II. Elements of Discretization 

The FVD schemes are derived from the integral form of a conservation law 

j (F • n) dT = JJ (f-S) dn, (1) 

r a 

where f is a forcing function independent of the solution, S' is a solution-dependent source function, Q is a control 
volume with boundary I’, n is the outward unit normal vector, and F is the flux vector. The general FVD approach 
requires partitioning the domain into a set of non-overlapping control volumes and numerically implementing equation 
(1) over each control volume. The focus of this paper is node-centered FVD schemes, where solution values are stored 
at the mesh nodes; however, the proposed analysis techniques can be straightforwardly applied to cell-centered FVD 
schemes and other discretization methods, e.g., finite difference or finite element. 

For 2D node-centered FVD schemes, the considered primal meshes are composed of triangular and quadrilateral 
cells; a median-dual partition can be constructed by connecting the mass centers of the primal-mesh cells with the 
midpoints of the surrounding edges (Figure 1). These non-overlapping control volumes cover the entire computational 
domain and compose a mesh that is dual to the primal mesh. 



Figure 1. Median-dual partition for node-centered finite-volume discretizations. Numbers at larger black circles denote grid nodes. Smaller 
black circles denote dual-mesh nodes. 

The discrete solution is represented as a piecewise polynomial function; the polynomials are defined either within 
primal or dual cells. The discretization is applied at a sequence of refined grids satisfying the consistent refinement 
property. The property requires the maximum distance across primal and dual cells to decrease consistently with 
increase of the total number of grid points, N. In particular, the maximum distance should tend to zero as TV” 1 / 2 in 2D 
computations, and as /V -1 / 3 in 3D computations. 

The main accuracy measure of any FVD scheme is the discretization error, Ed, defined as the difference between 
the exact continuous solution, Q, to the differential conservation law 

VF = / - 5 (2) 

and the exact discrete solution, Q h , of the discretized equations (1) 

E d = Q- Q h . (3) 

A common approach to evaluate the accuracy of a FVD scheme is to monitor the convergence of truncation errors. 
Truncation error, E t , measures the accuracy of the discrete approximation to the differential equations (2) and is 
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defined as the residual obtained after substituting the exact solution Q into the normalized discrete equations (1), 


Et M 


JJ (f h - S h (Q)) dS2- J (F h (Q) ■ n) dT 


where F h is a reconstruction of the flux F at the boundary 1’, Q is the measure of the control volume. 


(4) 


\n\ = JJ dn, ( 5 ) 

n 

f h and S h are, respectively, approximations of the forcing function f and the source function S on Cl, and the integrals 
are computed according to some quadrature formulas. 


III. Downscaling (DS) test 

Convergence of truncation errors has been applied as a FVD accuracy measure on both structured (regular) and 
unstructured (irregular) grids. 615 On structured grids, this approach is well justified, as the truncation error con- 
verges as 0(k p ) on sequences of consistently refined grids, where h is a characteristic meshsize and p is the design 
discretization-accuracy order of the method. However, the truncation-error convergence is often misleading for un- 
structured grids. Several studies 1115 noted that 2 T " , -ordcr convergence of truncation errors for some commonly used 
node-centered FVD schemes can be achieved only on grids with a certain degree of geometric regularity. Examples 
in 1 and in this paper show that the design-order discretization-error convergence can be achieved even when trunca- 
tion errors exhibit a lower-order convergence or, in some cases, do not converge at all. The DS tests and analytical 
estimates can be used to indicate whether the design order of the discretization accuracy can be achieved. 

The overall process is shown below, summarizing relations between error convergence orders predicted by analyt- 
ical estimates , observed in DS tests, and observed in global grid-refinement computations. The analytical estimates 
predict error convergence observed in DS tests; the estimates are conservative because they do not account for possible 
error cancellation occurring on regular grids. The DS test is designed to predict the convergence orders observed in 
global computations. The tests provide exact estimates for truncation error convergence, but do not account for possi- 
ble discretization-error accumulation. On truly unstructured multidimensional grids, both the analytical estimates and 
the DS-test estimates proved to be sharp. 


Analysis Estimates 

• conservative for discretization errors 

• conservative for truncation errors 

DS-Test Estimates 

• optimistic for discretization errors 

• exact for truncation errors 

Grid-Refinement Computations 


To apply a DS test, one first chooses a (manufactured) solution defined on the given computational domain. The 
associated forcing functions are found by substituting this solution into the continuous governing equations and bound- 
ary conditions. The DS test requires numerical computations on a sequence of contracted domains zooming toward a 
focal point within the original computational domain (Figure 2). The choice of the focal-point location can be varied to 
study a typical interior discretization, a boundary discretization, or a specific singularity. Because the number of points 
in the DS-test domains is held (approximately) fixed, one can study solutions on grids with characteristic meshsizes 
much smaller than those of practical computations with a globally-refined grid sequence. 
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(a) Scaled grid 



(b) Independent grid generation (accounting 
for curved physical boundary) 


Figure 2. Illustration of DS computational domains. Black bullets mark the focal points; white squares mark the interface between the 
interior and the DS-test domain. 


There are at least two possible strategies for grid generation on these contracted domains: (1) The first strategy 
is termed “scaled grid”; with this strategy, the first (coarsest) computational domain is defined as a subdomain of 
the investigated global mesh containing the focal point; other (finer) domains and their mesh patterns are derived by 
scaling down this first domain (e.g., repeatedly multiplying all the distances from the focal point by a given factor, 
say, 1/2 or 1/4). (2) An independent grid can be generated on each domain, assuming that the consistent refinement 
property is satisfied, i.e., the maximum distance across a grid cell is scaled down with the same rate as the diameter of 
the contracted domains. 

The scaled-grid approach (sketch (a) in Figure 2) is especially useful for studying interior discretizations and 
straight boundaries. It is impractical to implement it for studies near a general (discretely defined) curvilinear boundary 
because the physical boundary shape should be preserved on each grid in the DS sequence (sketch (c) in Figure 2). 
In the DS test, the studied FVD scheme is supplemented with a set of boundary conditions at the interface between 
the interior and the DS domain; overspecification from the known manufactured solution is a possible choice. The 
freedom to choose the focal point, the shape of domains, and the type of boundary conditions greatly simplifies DS 
testing. In studying accuracy of discrete boundary conditions, i.e., when the focal point is at the physical boundary, the 
physical conditions are implemented at the studied boundary surface; overspecification can be applied at the interface 
between the interior and the DS domain (see sketches in Figure 2). 

The DS test evaluates local truncation and discretization-error convergence orders by comparing errors obtained 
in computations on different scales. The DS convergence of truncation errors is a precise indication of the truncation- 
error convergence (in the Loo-norm) observed in global grid-refinement computations, as long as the DS test samples 
all representative regions. Global convergence in integral norms, e.g., L\ -norm, may be higher order because these 
norms are less sensitive to fluctuations occurring locally. 

On grids with an increased regularity, the DS test may overestimate convergence of the discretization errors. Some 
global error accumulation can occur on such grids; because of its local nature and overspecified interface boundary 
conditions, the DS test is incapable of accounting for the global error accumulation. Our experience shows that 
on truly unstructured multidimensional grids, errors produced locally dominate globally accumulated errors, and the 
discretization-error estimates obtained in a DS test become sharp. 
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IV. Windowing 


In large-scale computations, such as those compared at the Drag Prediction Workshops, 10 a collection of grids is 
generated for grid-refinement studies. In structured-grid settings, it is theoretically possible to generate a family of 
grids with similar properties such as grid topology, node distribution, cell regularity, etc.; in 3D practice, the generated 
unstructured grids rarely adhere to this similarity principle. In unstructured-grid settings, the underlying principles 
used in the generation are often heuristic and thus the concept of family-of-grids is harder to articulate. The only 
practical limitation to the analysis, as discussed subsequently, is that the consistent refinement property be satisfied for 
the grid sequences analyzed. 

In order to provide a framework for assessing practical large-scale grids, we introduce the concept of windowing. 
A window is an arbitrarily-shaped subdomain within the computational domain (see sketch in Figure 3) serving as a 
reference frame for either DS or grid-refined testing. The windows can encompass interior and boundary points; points 
of the computational domain outside of the window can be overspecified by the exact (manufactured) solution. The 
windows can be tailored to investigate the effect of interior discretizations, boundary discretizations, or singularities. 



Figure 3. Sketch of possible windows superimposed on an unstructured grid. Regions denoted by solid line are grid-scaled windows (with 
black focal points); regions denoted by dashed line are windows preserving body geometry (with gray focal points). 

For grid-refinement testing, the window is fixed and equations are defined on a sequence of refined grids generated 
within this window. In order to establish a meaningful convergence order, the consistent refinement property should 
be satisfied for all of the grids in all windows. The assessment of refinement consistency could be used during the grid 
generation process itself to provide more meaningful grids for convergence-order investigations. 

The DS tests correspond to solutions on a sequence of windows, each window scaled down in physical extent 
around the focal point. If the scaled windows include only interior points or lie along a (piecewise) straight boundary, 
the grids (of similar complexity) can be generated within the windows using scaling from a given mesh. For windows 
including points along a curved boundary, the boundary shape should be preserved at each scale, involving grid 
adjustments in each window. In any instance, the consistent refinement property should be satisfied implying the 
maximum distance across the grid cells is downscaled with the same rate as the window downscaling. 

When developing and maintaining simulation software, or when examining the suitability of an existing scheme to 
large-scale applications, monitoring of truncation-error order is an inexpensive way to verify the accuracy of a scheme. 
In such tests, the solutions need not be determined, only the residuals evaluated with the known (manufactured) 
solutions are examined. Even such simple computations may be tedious or misleading for large-scale grid-refinement 
studies; large amounts of data are involved and integral norms of truncation errors do not provide sufficient information 
to isolate the source of large errors and/or to determine if the large truncation error is detrimental for the discretization- 
error order properties. An attractive feature of DS windowing is that grids of very small sizes can be used to examine 
whether the grid solutions are within the asymptotic range of the order property. 

As regards the maintenance or development phases, the DS examinations can be performed for a range of situa- 
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tions, from simpler interior discretization tests for very regular grids to more complex boundary condition tests for 
very irregular or mixed-type grids. In particular, in the interior of a structured grid, the truncation errors predicted by 
the DS test and observed in grid-refinement tests are expected to satisfy the design-order property. Near the bound- 
aries, the expectations are lower, especially for the DS test and for the maximum of the truncation error observed in 
grid-refinement tests; in many situations, the truncation-error convergence order is allowed to be lower on a small 
subdomain without affecting the global discretization accuracy. The range of tests should, of course, also involve 
variations in solution parameters, such as Mach number or Reynolds number. As mentioned above, on truly unstruc- 
tured grids, the truncation error convergence is allowed to be 1 order lower than the discretization design order for 
inviscid-dominated equations and 2 orders lower, for viscous-dominated terms. Reynolds number controls the vari- 
ation between these two regions. In maintaining a simulation capability, it is envisioned that every equation set and 
every boundary condition should have a set of verification tests. 

Discretization error testing, either in grid-refinement or downscaling, is more expensive because discrete solutions 
must be computed. The DS tests are far simpler than calculations over the entire domain, because the number of 
points on each grid is kept (approximately) constant. The interpretation of the results is also simpler; in particular, if 
the discretization error convergence is slower than expected, this is an unambiguous indication of deficiencies in either 
formulation or implementation. 


V. Analysis of DS-test discretization and truncation errors 

The DS-tests analyses estimate the convergence orders of the discretization and truncation errors observed in DS 
tests. The analysis provides conservative estimates; typically on regular grids, the truncation-error convergence rates 
are higher because of local error cancellations. 

The discrete conservation-law equations are 

j> (F h (Q h ) • fi) (IT = jj (f h - S h (Q h )) (IQ, (6) 

r n 

where the flux reconstruction F h and the source approximation S h depend on the node -centered discrete solution, 
Q h ; outward normal n, the discrete force function f h , and all integrals are approximated with certain accuracy by 
predefined methods. Assuming the discretization error to be small comparing to the exact solution, Q, (\E ( j\ << |Q|), 
the discretization error can be evaluated as 

E d « J-\Q)R(Q), (7) 

where R (Q) and J ( Q ) are the residual and the Jacobian of the discrete conservation law, respectively; 


R(Q) = j| 

■'(Q) = oq 


- § (F h (Q) • n) dr + ff ( f h - S h (<?)) dn 


L r 
i 

ID 


§ (F h (0) • n) dr + JJ S h (Q) dn 
r n 


( 8 ) 


|r| = = ft d_1 is the measure of the control-volume boundary, where d is the space dimension and ft, is a 

characteristic diameter of the control volumes. Note, that in distinction from (4), R(Q) is the residual of an integral 
conservation law and as such is normalized with |F|, not |0|. The truncation error relates to R(Q) as 


E t = 


R(Q ) 

ft 


(9) 


A scalar (nonlinear) equation is called non-degenerated for a given manufactured solution, if the coefficients of 
the leading terms in the equation linearized around this solution do not all vanish. For FVD schemes discretizing 
non-degenerated scalar conservation laws with zero source function S , a simple rule for evaluating the asymptotic 
order, mj, of J~ 1 (Q) = 0(h mj ) is 


toj = (order of solution differentiation in F). (10) 

For general systems of nonlinear equations, estimating the asymptotic order of J~ 1 {Q) is a more complicated task. 
Analytically, the asymptotic order can be predicted by analyzing an equivalent linear operator , E(Q), that is derived 
from the J(Q) operator (8) by formally replacing all spatial derivative operators contributing to the Jacobian with 
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0(1 /h k ) terms, where k is the differentiation order. We will refer to these replacements as equivalent substitutions. 
After equivalent substitutions, E(Q) can be inverted to estimate the order of ■J~ 1 (Q). For basic fluid equations, this 
analytical evaluation can be performed with a relative ease. 1 For more complicated equations, the DS test can be 
expanded to evaluate J~ 1 (Q) numerically. 

The discretization-error order also depends on the R(Q ) approximation order, run. R(Q) = 0(h mR ). There 
are four main contributors determining the accuracy of R(Q): (1) flux reconstruction, (2) control-volume boundary 
approximation, (3) flux integration, and (4) source and forcing term integration. A p th - order accuracy for R(Q) 
(mu = p) can be achieved, if the following four sufficient accuracy conditions are satisfied. The conditions are 
formulated in terms of the maximal degree of a polynomial, for which the considered approximation provides the 
precise outcome. 

1. Flux reconstruction accuracy The flux F h should be reconstructed precisely at the control-volume 
boundaries for an analytical flux F represented by polynomials of the (p — l) th degree. 

2. Control-volume boundary accuracy The integral of the outward normal to the boundary should be 
computed precisely for each boundary segment represented as a (p — 1 ) th -order curve, i.e., its (piecewise) parametric 
formulation with respect to the curve length involves only polynomials of the (p — l) th degree. This contributor often 
plays an important role only near the physical boundaries of the computational domain, where the imposed geometry 
must be approximated. 

3. Flux integration accuracy The flux integration method should be exact for fluxes represented as polyno- 
mials of the (p — l) th degree. Note that this condition relates to the accuracy of the integration scheme, and assumes 
no reconstruction errors. One should distinguish between the global and local integration accuracies. The global 
p th - order integration accuracy is achieved when the numerical integration of polynomials of (p — l) th degree over the 
entire closed control-volume boundary is exact, i.e., provides the same result as the analytical integration. To achieve 
the local integration accuracy of p th order, the (p — l) th degree polynomials should be integrated exactly over each 
segment of the boundary. Local accuracy is more strict and implies global accuracy of the same order or higher; global 
p th - order flux integration accuracy is sufficient. 

4. Source and FORCING term integration accuracy The integral over the control volume should be com- 
puted precisely for the integrated function represented by a polynomial of the (p — 2) th degree. 

Each of these sufficient accuracy conditions can be evaluated separately and independently. These conditions are 
conservative (not necessary) because higher-order approximation to R ( Q ) can be achieved due to error cancellation, 
especially, on regular (structured) grids. For conservation laws, with J~ 1 (Q) = 0(h mj ) and R(Q) = 0(h mR ), the 
truncation and discretization errors are evaluated as 

E t = 0(h mR ~ 1 ) and E d = 0(h mR+mj ). (11) 


VI. Two-dimensional computations : Laplace equations 

As a model of the diffusion terms in the Navier-Stokes equations, we solve the two-dimensional Laplace equation, 

A U = f, (12) 

with Dirichlet boundary conditions, on a series of primal mixed-element grids composed of triangles and quadrangles. 
The exact solution and forcing term are taken as U = [(sin(7r;r)) 2 + (sin(7r?/)) 2 ]/2, / = — 27r 2 [l — (cos(7ra;)) 2 — 
(cos(7ry)) 2 ]. The solution is represented as a piecewise polynomial function, with polynomials defined at the primal 
cells; the conservation law 


t 


VU-ndT = 


fd.n 


(13) 


is enforced on node-centered control volumes constructed by the median-dual partition. 

With reference to Figure 4, the integral flux through the dual faces adjacent to the edge [ro, it] is computed as 


VU • ndr = VU" • n R + VU L • n L , 


( 14 ) 


ABC 


7 of 17 


American Institute of Aeronautics and Astronautics 




Figure 4. Illustration of reconstruction on mixed grids with median-dual partition. 


where riR and rtL are directed areas of the corresponding dual faces. The gradient is reconstructed separately at each 
dual face as follows. 

For the triangular element contribution, the gradient is determined from a Green-Gauss evaluation at the primal- 
grid element, 

VU l = VUoi 2 . (15) 

The gradient overbar denotes a gradient evaluated by the Green-Gauss formula on the primal cell identified by the 
point subscripts. With fully-triangular elements, the formulation is equivalent to a Galerkin finite element scheme 
with a linear basis function. 3 - 4 For the quadrilateral element contribution, the gradient is evaluated as 

VU fl = VU 03 4i+ Ul ~ U ° - VU0341 • e 01 eoi, ( 16 ) 

Lki- r ol J 

where r , is the coordinate vector of node r, and 


eoi 


ri - r 0 

hr - r o I 


(17) 


is the unit vector aligned with the edge [r 0 ,iT]. Note that for grids with dual faces perpendicular to the edges, the 
edge-gradient is the only contributor. This approach to the gradient reconstruction is used to decrease the scheme 
susceptibility to odd-even decoupling. 8 9 In all cases, the linear solution reconstruction leads to the l st -order flux 
(gradient) reconstruction accuracy. 



Figure 5. A typical mixed-element unstructured grid generated with random splitting and random perturbation of the underlying quadri- 
lateral grid. 

The asymptotic order of computed according to (10), is nij — 1. The volume integral of the forcing 

term in (13) is evaluated as the value at the node multiplied by the dual volume. Overall, R(Q) = 0(h). The analysis 
predicts the discretization error E r j - 0(h 2 ) and the truncation error E f = 0(1). 

Each grid is formed from an underlying structured quadrilateral grid (Figure 5). In terms of a polar, (r, 9), coordi- 
nate system, the grid extent is defined as 9 € [7t/6,7t/3] in the circumferential direction and r £ [1, 2.2] in the radial 
direction. The decision to split (or not to split) each structured quadrangle into triangles is determined randomly; 
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(a) DS test 


(b) Grid refinement computations 


Figure 6. Convergence of the discretization and truncation errors for the Laplace equation solved on irregular mixed-element unstructured 
grids. 


approximately half of the quadrilaterals are split. In addition, the interior grid points are perturbed from their original 
position by random shifts in the range [— \/2/6, \/2/6] of the local mesh size in the radial direction. The sequences 
of globally refined grids are generated with 2 n+3 + 1 points in both the radial and circumferential directions, where 
n = 0, 1, 2, 3, 4. The sequences of DS grids are generated from a grid with 17 points in both the nominal radial and 
circumferential directions and downscaled about the center of the domain by a factor 2~ a , where s = 0, 2, 4, 6, 8. 
The grid topology remains unchanged. 

The L \ norms of truncation and discretization errors are shown in Figure 6 versus an equivalent mesh size param- 
eter, taken as the L\ norm of the square root of the dual volume. Although not shown, error convergence rates in the 
Loo norm are the same as the Li-norm rates. The analysis predicts the DS-test error convergence rate precisely. In 
grid-refinement computations, the truncation errors remain 0(1) and the discretization errors converge with 2" d -order, 
precisely predicted by the DS test. The reason for the 0(1) convergence of truncation errors is grid irregularity stem- 
ming from the perturbation to the grid points and the usage of mixed grids. The literature frequently interprets 0(1) 
convergence of truncation errors on irregular grids as evidence of an inconsistent scheme; 6,1115 this example clearly 
shows that this is not the case. Note that the residual of the integral conservation law R(Q) is 0(h) in this example. 

VII. Two-dimensional computations : incompressible Euler equations 

In this section we consider incompressible inviscid equations. The source function S is assumed to be zero. 
Inviscid fluxes for conservation of mass and momentum are defined as 



/3u 


(3v 

F = fi + gj = 

u 2 + p 

i + 

uv 


uv 


v 2 +p 


where the vector of unknowns, Q = [u,v,p\, includes the Cartesian velocities and the pressure; (3 is an artificial 
compressibility parameter, 3 taken as /3 = 1 here. 

The median-dual partition is applied. At each dual control volume, a polynomial solution approximation is con- 
structed. The approximation is required to coincide with the discrete solution value at the central node ro; the poly- 
nomial coefficients are determined through a least-square procedure involving neighboring nodes, i.e., nodes linked 
by an edge to the central node r 0 . If the set of neighboring nodes is insufficient to determine, uniquely and stably, 
the polynomial coefficients, the set may be expanded to involve neighbors of neighbors. Three node-centered FVD 
schemes are considered: an edge -reconstruction scheme and two face-reconstruction schemes. 
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The asymptotic order of J 1 ( Q ) for the non-degenerated cases considered in this paper in m.j = 0 for all three 
equations. The effect of degeneration on convergence of the discretization errors is considered in previous work. 1 

A. Edge-reconstruction scheme 

This FVD scheme is commonly used in unstructured-grid computations 3 ' 4,s ' 9 - 1 1 The numerical upwind flux at the edge 
midpoint is approximated according to the Roe scheme, 3 ' 14 

(F h • ft) = i [(F 0 • ft) + (Fj • ft)] - ^ |A(Q)| (Qi - Q 0 ), (19) 

where, Qo and Qi are the “left” and “right” solution reconstructions at the edge midpoint derived from the linear ap- 
proximations at the control volumes centered at ro and ri, respectively (see Figure 4); Fo and Fi are the corresponding 
“left” and “right” unsplit fluxes; ft is the unit vector aligned with the combined directed area vector n = ni, + n R , 
where nL and n R are outward normal (directed-area) vectors of the left and right segments, respectively, with their 
amplitudes equal to the corresponding segment lengths; |A(Q)| is Roe’s approximate Riemann solver matrix com- 
puted for Q = |(Qi + Qo)- With a linear solution approximation at each control volume, this is a 2 order accurate 
flux reconstruction. 

The flux integration over the two segments of the control-volume boundary linked at the edge midpoint is approx- 
imated by multiplying F h computed at the edge midpoint with the amplitude of the combined-directed-area vector, 

| n | . According to the definition introduced in the sufficient flux-integration accuracy (condition 3 in Section V), this 
computationally efficient integration scheme possesses global 2 nti -order accuracy on triangular grids, 3 providing the 
exact solution for the conservation laws with linear fluxes. The scheme provides only l s4 -order local integration accu- 
racy. On general (unstructured) quadrilateral and mixed-element grids, the global flux integration accuracy deteriorates 
to I s4 order. 

The forcing term integration over the control volume is approximated as the node value multiplied by the volume 
Q | . This approximation is precise for the integral of a constant (0 th -order polynomial) forcing term. The overall 
approximation accuracy for R(Q) is R(Q) = 0(h 2 ) on triangular grids and R(Q) = 0(h ) on general quadrilateral 
and mixed-element grids, leading to estimates of Ed = 0(h 2 ) and E t = 0(h) for triangular grids and E r j = 0(h) 
and E t = 0(1) for general grids. 



Figure 7. Illustration for edge-reconstruction flux integration scheme near the boundary. 

The boundary is approximated as a piecewise straight line connecting the grid nodes located at the physical bound- 
ary. The control volume around a boundary node, ro, (see Figure 7) is closed with the boundary segments [r R , ro] and 
[ro, rc]. The straight-line approximation provides a 2 T " / -ordcr accuracy to the curved boundary segments connecting 
nodes ro and r 2 . However, the approximation provides only a l st -order accuracy for the ro7?F part of the bound- 
ary segment . The flux reconstruction and integration over the remaining segments are the same as for the interior 
edges (segments [r^,r R ] and [r R ,r c\ evaluate numerical fluxes at r R and ry - using the directed areas ri! and ri 4 , 
respectively). 

Boundary conditions are enforced weakly through the boundary fluxes. The inflow boundary conditions are im- 
plemented by evaluating velocity contributions to the boundary flux from prescribed conditions; the outflow boundary 
conditions are implemented analogously by evaluating pressure boundary-flux contributions from prescribed condi- 
tions. 
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On triangular grids, the integral flux through the boundary segment [rs,ro] is computed as (F h [ rs ro ] • 112), 
where ri2 is the directed area of [re, r 0 ], and F h - r;j ro , is evaluated at vp = 5/6ro + l/6r2. The solution components 
prescribed in the boundary conditions, e.g., zero normal velocity components at tangency, both velocities at inflow, 
or the pressure at outflow, are specified at vp from the known exact solution; other components are interpolated to 
r p from the endpoints of the segment [r2,ro]. The integral flux through [r (l , re] is computed analogously. In spite 
of providing only l st -order local accuracy for integrated fluxes through [r^, r 0 ] and [r 0 , rc], the coefficients 5/6 and 
1 /6 lead to cancellation of l st -order errors, providing zero residuals for conservation laws with linear fluxes and, thus, 
supporting the global 2"' , -ordcr flux integration accuracy. 4 

On mixed-element and irregular-quadrilateral grids, the 2" d -order flux integration accuracy is not recovered. For 
triangular cells adjacent to the boundary, the (5/6, 1/6) rule is still enforced; on the boundary segments of quadrilateral 
cells, the fluxes are evaluated at the segment midpoint, r p = 3/4r 0 + l/4r 2 . The overall boundary approximation 
accuracy for R(Q) is consistent with the accuracy of the interior FVD scheme, 2 nd order for triangular grids and 1 st 
order for other irregular grids. 

Reliance on error cancellation in providing design-order discretization accuracy may lead to some counterintuitive 
phenomena. In particular, improving the accuracy of boundary conditions in a way that is not compatible with the flux 
computations through interior segments can upset the error cancellation balance and, in fact, worsen approximation 
accuracy of R(Q), at least, locally. In some cases, this local imbalance is compensated by contributions from other 
boundary segments. The imbalance is especially prominent at the corners where the boundary segments represent 
different types (e.g., tangency and inflow) of boundaries. 

B. Face-reconstruction schemes 

In this section, we describe two face-reconstruction FVD schemes that employ the median-dual partition and provide 
2 nd and 3 rd order accuracies on general mixed grids. Similar to the edge-reconstruction scheme, the linear and 
quadratic polynomials are defined at dual control volumes and coincide with the solutions at the grid nodes. The 
polynomial coefficients are defined in a least-square procedure involving solutions at the neighboring nodes. In the 
interior, piecewise straight dual boundaries imply a precise representation for n. The schemes described are quite 
similar to the scheme described by Delanaye and Liu 7 for cell-centered discretizations. As noted there and also shown 
in the appendix of the paper introducing DS, 1 the operation counts significantly favor the cell-centered approach for 
flux integration in three dimensions. 

The distinguishing feature of these face-reconstruction schemes is the flux integration procedures ensuring the 
designed local integration accuracy. Extension of 2 m, -ordcr accuracy to mixed grids can be achieved with linear poly- 
nomials, a modified straight-segment tangency boundary approximation (see Figure 10 and discussion in Section C), 
and a flux integration scheme providing local 2” d -order accuracy. Accuracy of 3 rd order requires quadratic polynomi- 
als for flux reconstruction, a quadratic fit to the curved tangency boundary, an integration scheme with local 3 rd -order 
accuracy, and a 2 nd -order accurate scheme for integrating source and force terms. 

In the implemented 2D version, a polynomial flux is defined at each segment of the dual control-volume boundary. 
In (19), the “left” and “right” variables, Qo and Qi, are represented by the basic polynomials defined at the adjacent 
control volumes. The dissipation matrix, |A(Q)|, is computed for Q, defined as the average of Q 0 and Q t evaluated 
at the segment midpoint. The scheme with linear polynomials provides 2 nd -order accurate flux reconstruction; fluxes 
reconstructed with quadratic polynomials are 3 rd -order accurate. 

The inner product of the polynomial flux vector and the outward unit normal vector is integrated over the segment 
with a numerical Gauss-Legendre quadrature formula employing 1 point for the linear face-reconstruction scheme, 
and 2 points for the quadratic face-reconstruction scheme. As a convenient debugging tool, one can add one Gaussian 
point to have total of 2 and 3 points, respectively, per integration segment. With these improved accuracy integrations, 
the FVD method should provide zero residuals for linear and quadratic manufactured solutions. 

The forcing term integration method for the linear face-reconstruction scheme is the same as that for the edge- 
reconstruction scheme. For the quadratic face-reconstruction scheme, a linear (or piecewise linear) approximation to 
the forcing term is constructed and integrated at each control volume providing a 3 rd -order accurate approximation. 
The overall approximation accuracy for R(Q) is 2 nd order for linear polynomials and 3 rd order for quadratic poly- 
nomials and does not deteriorate for general unstructured grids. The predictions for the error convergence rates are 
Ed = Oih 2 ), E t = 0(h) for the linear polynomials and Ed = 0(h 3 ), E t = 0(h 2 ) for the quadratic polynomials. 

The boundary conditions for the face-reconstruction discretizations are enforced weakly, through boundary fluxes. 
At inflow and outflow, the computational domain is bounded by a piecewise straight boundary; the solution compo- 
nents are either specified from the exact solution (the velocity is specified at inflow and the pressure is specified at 
outflow) or represented by the polynomial approximation defined at the adjacent control volume. At tangency, the 
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2 -order accuracy can still be achieved with enforcing tangency over piecewise linear segments. (See discussion in 
Section C.) For the 3 rY/ - order accuracy, the physical boundary should be approximated quadratic ally. In a general case, 
where the analytical shape of the boundary is unknown, the boundary should be represented as a piecewise polynomial 
curve providing the required accuracy for the boundary shape. 


C. Numerical tests for non-degenerated flows 

Numerical tests presented in this section are performed for 2D inviscid incompressible non-lifting flow around a 
cylinder of unit radius centered at the origin. The flow is described by the conservative equations (1) with zero source 
and forcing terms and fluxes defined in (18). The analytical solution for this problem is well known. 1 



(a) Random triangular 

Figure 8. Typical triangular and mixed-element 
the cylinder. 




(b) Random quadrilateral (c) Random mixed 

unstructured grids defined on a computational domain shifted away from the surface of 


Interior inflow/outflow domains The first set of tests is performed on a computational domain shifted away 
from the surface of the cylinder: 1.5 < r < 4, 27r/3 < 9 < 47t/ 3. Seven formulations are studied: the edge- 
reconstruction FVD scheme on random triangular, random quadrilateral, and random mixed-element grids; and face- 
reconstruction FVD schemes on random quadrilateral and random mixed-element grids. Examples of unstructured 
grids derived from an underlying structured grid are shown in Figure 8. For triangular and mixed-element grids, 
randomization is introduced through random splitting (or not splitting) of structured quadrilateral cells. Each cell has 
equal probabilities to introduce either of the two diagonal choices or, for mixed-element grids, no diagonals. For 
random quadrilateral grids, structured interior nodes are perturbed from their original position by random shifts in the 
range [— \/2/6, v/2/6] of the local mesh size in both the radial and circumferential directions. 

For each formulation, grid refinement and DS tests are performed. In grid refinements, the underlying structured 
grid is refined by doubling the number of intervals in the radial and angular directions. In the DS test, the coarsest 9x9 
grid is scaled down around the point r = 2.75, 0 = tt by multiplying all angular and radial differences from this point 
by a factor of 0.5. Randomization is introduced independently on each scale. The inflow boundary conditions are 
enforced at the boundary corresponding to the external radius; outflow conditions are enforced at all other boundaries. 
Table 1 summarizes the convergence of discretization and truncation errors observed in these tests. The convergence 
orders are the same in all norms and for all variables and equations. The results confirm analysis predictions and 
capabilities of the DS test to provide sharp estimates for error convergence in grid refinement computations. The 
convergence rates observed for the edge-reconstruction scheme on random triangular and random quadrilateral grids 
are consistent with previous results. 2 Although not shown, we have implemented a central version of the edge- 
reconstruction scheme. The observed convergence rates were identical with the rates shown in Table 1 . The results 
contradict Svard and Nordstrom, 1 1 where 0(1) convergence of discretization errors on random quadrilateral grids with 
a central-difference scheme was reported. 

Figure 9 shows L \ -norms of discretization error in it during grid refinement using the random quadrilateral grids 
for the three schemes; the behavior in u is typical of all variables. The edge-reconstruction scheme is T s, -ordcr on the 
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Formulation 

DS 

test 

Grid-refinement 

computations 

Discr. Error 

Trunc. Error 

Discr. Error 

Trunc. Error 

Edge-reconstruction, 
random triangular grid 

0(h 2 ) 

0(h) 

0(h 2 ) 

0(h) 

Edge-reconstruction, 
random quadrilateral grid 

0(h) 

0(1) 

o(h ) 

0(1) 

Edge-reconstruction, 
mixed-element grid 

0(h) 

0(1) 

o(h ) 

0(1) 

Linear 

face-reconstruction, 
random quadrilateral grid 

o(h 2 ) 

0{h) 

o(h 2 ) 

0(h) 

Linear 

face-reconstruction, 
mixed-element grid 

o(h 2 ) 

0(h) 

o(h 2 ) 

o(h) 

Quadratic 

face-reconstruction, 
random quadrilateral grid 

o(h 3 ) 

o(h 2 ) 

o(h 3 ) 

o(h 2 ) 

Quadratic 

face-reconstruction, 
mixed-element grid 

o(h 3 ) 

o(h 2 ) 

o(h 3 ) 

o(h 2 ) 


Table 1. Convergence of discretization and truncation errors for various unstructured grid formulation of the 2D inviscid incompressible 
equations on an inflow/outflow computational domain. 


9x9 grid 



Figure 9. Li -norms of grid-refined discretization errors in u using the random quadrilateral grids. 


finer grids. The linear and quadratic face-reconstructions recover design order. 
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TANGENCY BOUNDARY For the edge-reconstruction FVD scheme on triangular grids, local accuracy deterioration 
occurs if a curved tangency boundary is approximated by straight segments linking primal-mesh nodes located at the 
physical boundary. Sketch (a) of Figure 10 illustrates this approximation: the straight segments [r 2 ,ro] and [ro,r< 4 ] 
approximate a curved boundary, the points and r q are the segments’ midpoints, and the arrow indicates the local 

flow velocity. For the edge-reconstruction scheme, discrete tangency is enforced over the straight segments [ro, r/ 3 ] 
and [ro, r c\- The exact continuous solution satisfies the tangency condition at the actual curved boundary, not at the 
straight boundary segments. As estimated analytically 1 and confirmed numerically in Section C, this approximation is 
l st -order for each segment of the tangency boundary. 



(a) Basic median-dual approximation 


(b) Modified approximation 


Figure 10. Straight-segment approximations to curved tangency boundary. 



x 


Figure 11. Random triangular grid around the top of the cylinder. 

To illustrate the effect of a straight-segment approximation to a curved tangency boundary, a sequence of random 
triangular grids is generated at the top of the cylinder (1 < r < 2.2, tt/3 < 9 < 2tt/ 8) and used in computations 
with the edge-reconstruction FVD scheme; a grid example is shown in Figure 11. Figure 12 illustrates convergence 
of the L i norm of truncation and discretization errors in grid-refinement and DS tests. Two DS tests are performed, 
each with the focal point at the top surface of the cylinder. The first uses overspecification at all boundaries except 
the interior tangency nodes. The second replaces the overspecification along one boundary with the physical inflow 
boundary condition (see sketches in Figure 12). 

For the interior-tangency nodes, there are two tangency segments; the errors at these segments contribute to R(Q) 
with opposite signs and, at least partially, compensate each other. Because of this compensation, the accuracy deteri- 
oration does not affect interior-tangency nodes on the grids with (nearly) uniform boundary node distributions. The 
2 nd -order convergence of discretization errors and the L** -order convergence of truncation errors demonstrated in the 
interior-tangency DS test confirm this conclusion. However, in the corners or/and at the interfaces between clusters 
of boundary nodes with different topology/spacing, the compensation does not occur. The accuracy deterioration is 
clearly observed in the DS test performed with the inflow/tangency boundary conditions. The convergence of the L\ 
norms of the errors in the grid-refinement test is the same as convergence shown in the interior-tangency DS test and 
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(a) DS test: interior tangency boundary conditions 


(b) DS test: inflow/tangency boundary conditions 
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Truncation Error : x- momentum 
Truncation Error : y- momentum 
Discretization Error : p 
Discretization Error : u 
Discretization Error : v 


(c) Grid refinement 


Figure 12. Convergence of the L\ norm of truncation and discretization errors observed in DS and grid-refinement tests for the edge- 
reconstruction FVD scheme. The tests are performed on random triangular grids surrounding the top tangency boundary of the unit 
cylinder. The open square symbols in the sketches denote overspecification for the DS tests. 


is not affected by local accuracy deterioration in the corners; although not shown, the grid-refinement convergence of 
the Lqo norm of errors is similar to the convergence demonstrated in the DS test with the inflow/tangency corner. 

The linear face-reconstruction scheme possesses the flexibility to recover the 2" d -order accuracy with a straight- 
segment tangency boundary approximation. The required modifications is illustrated in the sketch (b) of Figure 10: 
the dual r b and r c are moved to the boundary. Figure 13 shows convergence of the L-\ norm of truncation and dis- 
cretization errors observed in DS tests for the linear face-reconstruction FVD scheme with the two types of straight- 
segment approximation to the curved tangency boundary: the basic median-dual approximation and the modified 
approximation. The tests are performed in the setting similar to the inflow/tangency boundary condition used with the 
edge -reconstruction scheme. The results confirm that the modified approximation, 2 rad -order convergence is achieved. 
Although not shown, with the quadratic face-reconstruction scheme and a quadratic fit to the curved tangency bound- 
ary, discretization errors converge with the 3 rd order. 

VIII. Three-dimensional computations : Incompressible turbulent equations 

In the full paper, we will show both DS and grid-refinement tests for the incompressible turbulent flow equations 
on unstructured grids. The tests show that the convective operator in the turbulence equation can be evaluated with 
the first-order accuracy yet second-order discretization accuracy in the meanflow solutions is recovered. This is a 
rather unexpected result, but it is easily seen to be true based on the Laplace equation example given above as an 
example of the diffusion terms in the Navier-Stokes equations. The upshot of solving the turbulence equation is to 
produce an eddy viscosity term, which is simply a coefficient in the flux gradient term. As shown above, the l st -order 
accurate flux-gradient approximation is sufficient for solutions with 2 -order discretization accuracy. The results 
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Figure 13. Convergence of the L\ norm of truncation and discretization errors observed in DS tests for the linear face-reconstruction FVD 
scheme with two types of straight-segment approximation to the curved tangency boundary. 


are necessarily a function of the Reynolds number. With a turbulence and meanflow variation prescribed through the 
method of manufactured solutions, very low Reynolds number (turbulent diffusion-dominated) flows show second 
order for both the turbulent and the meanflow equations. Higher Reynolds number flows, where convection is more 
important, show first-order discretization accuracy in the turbulence but second-order discretization accuracy in the 
meanflow solutions. 


IX. Concluding Remarks 

New concepts to verify large-scale computational methods using unstructured grids are presented. The method 
introduces computational windows, easily constructed within a collection of grids or a single grid, to examine or- 
der properties of the scheme over the entire domain or in regions isolating the interior discretization, the boundary 
discretization, or singularities. A major component of the methodology is the downscaling (DS) test, introduced ear- 
lier for studying the convergence rates of truncation and discretization errors of finite-volume discretization (FVD) 
schemes on general unstructured grids. Demonstrations of the method are shown, including a comparative assessment 
of the accuracy from a common edge-reconstruction FVD scheme, known to be second-order accurate for inviscid 
equations on triangular grids, with alternative second and third order accurate face-reconstruction schemes on gen- 
eral mixed grids. Improvements to the local accuracy at intersections of tangency and inflow/outflow boundaries are 
identified and repaired using DS tests. Additionally, analysis and computations are presented demonstrating that the 
usual discretization of first-order convection in the turbulence equations does not degrade the order of accuracy of the 
meanflow equations. 
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